Swing‐out opening of stromal interaction molecule 1

Abstract Stromal interaction molecule 1 (STIM1) resides in the endoplasmic reticulum (ER) membrane and senses luminal calcium (Ca2+) concentration. STIM1 activation involves a large‐scale conformational transition that exposes a STIM1 domain termed “CAD/SOAR”, ‐ which is required for activation of the calcium channel Orai. Under resting cell conditions, STIM1 assumes a quiescent state where CAD/SOAR is suspended in an intramolecular clamp formed by the coiled‐coil 1 domain (CC1) and CAD/SOAR. Here, we present a structural model of the cytosolic part of the STIM1 resting state using molecular docking simulations that take into account previously reported interaction sites between the CC1α1 and CAD/SOAR domains. We corroborate and refine previously reported interdomain coiled‐coil contacts. Based on our model, we provide a detailed analysis of the CC1‐CAD/SOAR binding interface using molecular dynamics simulations. We find a very similar binding interface for a proposed domain‐swapped configuration of STIM1, where the CAD/SOAR domain of one monomer interacts with the CC1α1 domain of another monomer of STIM1. The rich structural and dynamical information obtained from our simulations reveals novel interaction sites such as M244, I409, or E370, which are crucial for STIM1 quiescent state stability. We tested our predictions by electrophysiological and Förster resonance energy transfer experiments on corresponding single‐point mutants. These experiments provide compelling support for the structural model of the STIM1 quiescent state reported here. Based on transitions observed in enhanced‐sampling simulations paired with an analysis of the quiescent STIM1 conformational dynamics, our work offers a first atomistic model for CC1α1‐CAD/SOAR detachment.

: CC1α1-CAD/SOAR interface residues sorted by mean reweighted contact frequency ⟨ω ij ⟩. ⟨·⟩ refers to the average taken over 12 independent metadynamics runs. WT (5) K413S (6) E255S (6) K377S (5) A B D E F C Figure S3: (A) Zoom-in on crucial interaction partners forming CC1α1-CAD/SOAR salt bridges. (B) Mean electrostatic interaction energy for the 10 strongest CC1α1-CAD/SOAR salt bridges in the STIM1 monomer and dimer, respectively. For the monomer, interaction energies were calculated over 500 ns of conventional MD. For the dimer, mean interaction energies were calculated for three 500 ns replicas. Here, error bars denote the standard deviation over the individual replicas. (C) Number of hydrophobic residues near E255, E370, K377 and K413, and the average number of hydrophobic neighbours for all charged residues in our model. Note the above-average number of hydrophobic residues surrounding K377 and K413, which suggest that mutations K377A and K413A effectively trade electrostatic CC1α1-CAD/SOAR interactions for hydrophobic ones, overall preserving STIM1 function. The number of hydrophobic neighbours was calculated using a distance cutoff of 9Å. Distances were calculated between the Cβ atom of the charged residue and all heavy sidechain atoms of hydrophobic residues. (D-F) Patch clamp time series of whole-cell inward currents for the STIM1 WT together with ten different mutants, H259A/S, E370A/K, K377A/S K413A/S and E255A/S. Positions 377 and 413, which are surrounded by a large number of nearby hydrophobic residues (C), induce constitutive STIM1 activation when mutated to serine, but preserve storeoperated STIM1 function when mutated to alanine.   Figure S6: Overview of CRAC channel hallmarks, i.e. activation time and maximum current density I max , for all mutants tested in patch clamp experiments. Mutants that preserve store-operated STIM1 function are marked in green. Note that in line with their inability to strongly affect STIM1 function in the patch clamp experiment, all the corresponding mutated positions exhibit a low simulated binding score S i ranging from 0.01 to 0.62. Mutants marked with a "+" are constitutively active at time t = 0. Figure S7: Pairwise RMSD between structures from the bound state ensemble (defined as G < 2 kcal/mol). Overall, constitutively active STIM1 mutants (orange) show a broadened bound-state free energy minimum ( Figure S4). For these mutants, pairwise RMSD calculated for all bound-state conformations is on average higher than for the WT, i.e. a more diverse ensemble of conformations is accessible in the bound state. To faithfully reflect structural diversity of different CC1α1-CAD/SOAR bound states, RMSD was calculated based on Cα positions of interface residues 234-270 and 343-440.

M244S
L251S K413S L416S E370A I409S Figure S8: For selected mutants, the difference in binding score S i with respect to the WT is color-coded onto each binding residue. Red (blue) colors indicate stronger (weaker) binding compared to the WT. For each mutant, the mutated position is highlighted using space-filling representation (arrows).    Table S2: Residue pairs used for Figure 5B. Residues in the second STIM1 monomer are denoted by a prime. smFRET distances were taken from [1]. As in [1], 10Å were subtracted from each smFRET-derived distance to account for fluorophore linkers.
Simulated data refers to center of mass distances averaged over three independent MD replicas, each with a length of 500 ns. Figure S13: For each residue pair with available smFRET distance data, the difference ∆d between our modelled distance d MD and the experimentally determined distance d exp is shown (time evolution is color-coded from blue to yellow). Initial distances d MD (t = 0) obtained from our docked model, before simulating molecular dynamics, are denoted by red circles. Modelled distance ∆d(t) is smoothed by a time average over all time steps t i < t, i.e. ∆d(t) = mean out of 36 residue pairs, simulated distances relax towards d exp over the course of the simulation. Note that those residue pairs where the deviation between simulated and experimental distances increases over the course of the simulation, 309:309' and 312:312', are situated at the tip of the CC1α1,2 hairpin, which is an especially mobile region (see Figure S15).    In patch clamp experiments, one probes the free energy difference ∆G = G act -G quies between the activated and the quiescent STIM1 states, which can be dissected into three contributions, where ∆G unb +∆G unb→open is the free energy of STIM1 opening illustrated in Figure  6D and ∆G open→act describes the activation of Orai1. From our simulations we extract ∆G unb rather than ∆G. We assume that mutations in the CC1α1-CAD/SOAR binding interface only affect the free energy of unbinding, ∆G unb , i.e. ∆G unb→open and ∆G open→act are constant for different mutations (see Figure 6D). The overall equilibrium constant of STIM1 activation may be decomposed according to where the equilibrium constant of STIM1 opening/elongation reads and that of STIM1 activation by coupling to Orai1 is given as The experimentally determined CRAC channel current before store depletion, I(t = 0), is proportional to the concentration of bound STIM1-Orai1 complexes, where n denotes the cooperativity parameter [3][4][5], which takes into account that more than one STIM1 dimer binds to Orai1 and that successive binding events influence each other [6,7]. By we obtain where c 0 is a proportionality constant that also accounts for K STIM1 act , [STIM1 total ] and [Orai1] and c 1 = exp ∆G unb→open RT .

Model construction
For docking CC1α1 (PDB id 6YEL [10]) to CAD/SOAR using the Haddock 2.4 webserver [11,12], unambiguous restraints were added to reinforce pairing between BS 3 crosslinked residues: distances between residue pairs 413:243, 413:246, 386:238, 384:238 and 246:384 were restrained to values below 10Å with default energy constant scaling for the four Haddock stages. Clustering was performed based on the fraction of native contacts. All other docking parameters were kept at their default values. To select the starting model for our MD simulations, clusters were sorted by the Haddock score of the best-scoring structure within each cluster (instead of their average score; see Figure S18A).
For the assignment of protonation states, we calculated electrostatic potentials using Tapbs [13]. Protonation states were assigned according to protonation patterns sampled using Karlsberg 2.0 [14]. All titratable residues were found to occupy their standard protonation states, except for histidines, which showed varying patterns of protonated or neutral tautomeric states. The water solvating the final monomeric STIM1 model box was 160 × 110 × 95Å 3 in size, containing the protein, 51186 water molecules and 8 chloride counterions. To suppress rotation of our model, the positions of backbone atoms in CC1α1 were restrained during equilibration.
For the calculation of electrostatic interaction energies for CC1α1-CAD/SOAR salt bridges, we used gRINN [15] with the same settings employed in Namd [16] non-bonded interaction calculation.

Metadynamics methods
The number of contacts was calculated with the contacts collective variable implemented in the colvars module [17], using parameters cutoff = 6Å, expNumer = 6 and expDenom = 12. During out metadynamics runs, CC1α1 backbone positions were restrained to provide a firm basis from which CAD/SOAR is detached. All simulations were propagated up to a specific cutoff time, which was determined by calculating the number of contacts, N c , between CC1α1 and CAD/SOAR and terminating the run once N c fully zeroes out for at least 200 ps. This procedure ensures that all simulations stop just when the initial CC1α1-CAD/SOAR binding minimum is fully filled up by the bias potential. Secondly, since the contact frequency relies on the notion of a total run length t c (see below), this cutoff scheme facilitates a consistent definition of t c across all runs. Two-dimensional free energy surfaces were obtained by integrating out the contact collective variable. The CC1α1-CAD/ SOAR binding free energy is then given by the "depth" of the free energy well, i.e. G min -G d∞ , where d ∞ is a large, unsampled CC1α1-CAD/SOAR distance where the calculated free energy profile is flat.
Reweighting of non-biased observables was done using the balanced exponential reweighting scheme by Schäfer and Settanni [18], with weights where s(r) denotes a collective variable at point r in coordinate space, V (s(r), t) is the time-dependent metadynamics bias potential, β = 1 k B T and ⟨·⟩ denotes the average over the collective variable s.
The biased contact frequency ω b ij was calculated directly from our metadynamics runs as for a pair of binding residues i and j, one of which is in CC1α1 and one in CAD/SOAR. Using the weights w bex t , the reweighted contact frequency is given by Here, ∆t is the time interval at which atomic coordinates were stored (20 ps), θ denotes the step function with cutoff distance d c = 5Å and d ij (t) is the lowest interatomic distance between residues i and j at frame t. were chosen as active residues. On grounds of Figure 5-figure supplement 2 and Figure 2B in reference [1], position 371 was added as a further active residue and the distance between residues 400 and 400' was restrained to values below 50Å. To stabilize parallel orientation of the two monomers, the distance between residues 218 and 218' was restrained to values below 90Å. All other docking parameters were set as described in section 4.1. Comparing with smFRET-derived distances reported in [1] showed that the angle between our two monomers needed to be corrected. To this end, we targeted and subsequently restrained distances between residue pairs 309:309', 312:312', 388:388', 389:389' and 399:399' using harmonic restraints with target distances of 35Å, 37Å, 48Å, 46Å, and 42Å (derived from peak FRET values shown in Figure 2B in [1] and using Supplementary file 1 in [1]) while treating helical elements as rigid bodies by restraining dihedral angles [20]. From the 10 ns restrained MD run, we picked the frame with the lowest inter-monomeric interaction energy (calculated using Charmm [21]). Using this structure, we assigned protonation states using Tapbs and Karlsberg 2.0 as described above, but assigning a dielectric constant of 2 to a slab of 27Å thickness to model the membrane as a separate dielectric medium. All titratable residues were found to occupy their standard protonation states, except for histidines, which showed varying patterns of protonated or neutral tautomeric states. The anisotropic network model was constructed with ProDy [22] with a springs connecting all atoms within a 15Å radius and a unit spring constant.

FRET Microscopy
The apparent FRET efficiency E app was calculated after threshold determination and background signal subtraction. This was done on a pixel-to-pixel basis using a custom program [23] integrated into MATLAB (v7.11.0, The MathWorks, Inc., MA, USA) [24] and implementing the method described in [25] with a microscope-specific constant G parameter of 2.75. The Pearson correlation coefficient (R value) was used to measure the strength of the linear association/co-localization between STIM1 and Orai1 variants, where a value of R = 1 signifies perfect positive correlation/colocalization.

Electrophysiology
For store-dependent currents, the initial current amplitude recorded during the voltage ramp applied immediately after achieving the whole-cell configuration was subtracted from all subsequent current amplitudes. For constitutively active currents, a current amplitude obtained during a voltage ramp applied after perfusion with 10 µM La 3+ at the end of the experiment was used instead. Individual experiments were normalized by dividing all current amplitudes by the whole-cell capacitance, resulting in current density in picoampere per picofarad (pA/pF), to compensate for potential cell dimension-related effects and make individual experiments comparable. Individual time-series were aligned such that they reach maximum activation at the same time.